{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:58:06.468311",
     "start_time": "2016-08-17T07:58:04.883583"
    },
    "autoscroll": "json-false",
    "collapsed": false,
    "ein.tags": [
     "worksheet-0"
    ]
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import sys \n",
    "from os import getcwd, path\n",
    "sys.path.append(path.dirname(getcwd()))\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sb\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pprint\n",
    "from cohorts.functions import *\n",
    "import lifelines as ll\n",
    "import patsy\n",
    "import functools\n",
    "import survivalstan\n",
    "from cohorts.utils import strip_column_names\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## set seeds for stan & rngs, to aid in reproducibility\n",
    "## (note: only reproducible within a machine, not across machines)\n",
    "\n",
    "seed = 91038753\n",
    "import random\n",
    "random.seed(seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:58:06.519332",
     "start_time": "2016-08-17T07:58:06.499244"
    },
    "autoscroll": "json-false",
    "collapsed": false,
    "ein.tags": [
     "worksheet-0"
    ]
   },
   "outputs": [],
   "source": [
    "from utils import data\n",
    "from utils import paper\n",
    "from utils.extra_functions import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# prep data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:58:11.840853",
     "start_time": "2016-08-17T07:58:06.520959"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inner join with tcr_peripheral_a: 25 to 25 rows\n",
      "inner join with ensembl_coverage: 25 to 25 rows\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inner join with tcr_peripheral_a: 25 to 25 rows\n",
      "inner join with ensembl_coverage: 25 to 25 rows\n",
      "{'dataframe_hash': 7698303973572390439,\n",
      " 'provenance_file_summary': {u'cohorts': u'0.4.0+3.gda968fb',\n",
      "                             u'isovar': u'0.0.6',\n",
      "                             u'mhctools': u'0.3.0',\n",
      "                             u'numpy': u'1.11.1',\n",
      "                             u'pandas': u'0.18.1',\n",
      "                             u'pyensembl': u'1.0.3',\n",
      "                             u'scipy': u'0.18.1',\n",
      "                             u'topiary': u'0.1.0',\n",
      "                             u'varcode': u'0.5.10'}}\n"
     ]
    }
   ],
   "source": [
    "cohort = data.init_cohort(join_with=[\"ensembl_coverage\",\"tcr_peripheral_a\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:21.487992",
     "start_time": "2016-08-17T07:58:11.931783"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inner join with tcr_peripheral_a: 25 to 25 rows\n",
      "inner join with ensembl_coverage: 25 to 25 rows\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inner join with ensembl_coverage: 25 to 25 rows\n"
     ]
    }
   ],
   "source": [
    "def tcell_fraction(row):\n",
    "    return row[\"T-cell fraction\"]\n",
    "\n",
    "def peripheral_clonality_a(row):\n",
    "    return row['Clonality']\n",
    "\n",
    "cols, d = cohort.as_dataframe([snv_count,\n",
    "                               missense_snv_count,\n",
    "                               neoantigen_count,\n",
    "                               expressed_exonic_snv_count,\n",
    "                               expressed_missense_snv_count,\n",
    "                               expressed_neoantigen_count,\n",
    "                               exonic_snv_count,\n",
    "                               peripheral_clonality_a,\n",
    "                               tcell_fraction,\n",
    "                               ],\n",
    "                              rename_cols=True,\n",
    "                              return_cols=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:21.537104",
     "start_time": "2016-08-17T07:59:21.513702"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['snv_count',\n",
       " 'missense_snv_count',\n",
       " 'neoantigen_count',\n",
       " 'expressed_exonic_snv_count',\n",
       " 'expressed_missense_snv_count',\n",
       " 'expressed_neoantigen_count',\n",
       " 'exonic_snv_count',\n",
       " 'peripheral_clonality_a',\n",
       " 'tcell_fraction']"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cols"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## construct/rescale variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:21.559557",
     "start_time": "2016-08-17T07:59:21.539192"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## add/modify count variables\n",
    "d['nonexonic_snv_count'] = d.snv_count - d.exonic_snv_count\n",
    "cols.append('nonexonic_snv_count')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:22.132709",
     "start_time": "2016-08-17T07:59:22.081884"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## create 'observed', log-transformed & centered versions of variables (not normalized by MB)\n",
    "for col in cols:\n",
    "    observed_col = 'observed_{}'.format(col)\n",
    "    log_col = 'log_{}'.format(col)\n",
    "    log_col_centered = 'log_{}_centered'.format(col)\n",
    "    log_col_rescaled = 'log_{}_rescaled'.format(col)\n",
    "    d[observed_col] = d[col]*d['mb']\n",
    "    d[log_col] = np.log1p(d[observed_col])\n",
    "    d[log_col_centered] = d[log_col] - np.mean(d[log_col])\n",
    "    d[log_col_rescaled] = d[log_col_centered]/np.std(d[log_col_centered])\n",
    "\n",
    "## save key vars in list for future use\n",
    "vars_centered = ['log_{}_centered'.format(col) for col in cols]\n",
    "vars_rescaled = ['log_{}_rescaled'.format(col) for col in cols]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:22.174374",
     "start_time": "2016-08-17T07:59:22.134744"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## construct new variables for key ratios / comparisons\n",
    "\n",
    "# what proportion of X are expressed?\n",
    "d['exonic_expression_ratio'] = d.expressed_exonic_snv_count / d.exonic_snv_count\n",
    "d['missense_expression_ratio'] = d.expressed_missense_snv_count / d.missense_snv_count\n",
    "d['neoantigen_expression_ratio'] = d.expressed_neoantigen_count / d.neoantigen_count\n",
    "\n",
    "# d['expressed_missense2snv_ratio'] = d.expressed_missense_snv_count / d.expressed_snv_count\n",
    "d['missense2exonic_snv_ratio'] = d.missense_snv_count / d.exonic_snv_count\n",
    "d['expressed_neoantigen2missense_ratio'] = d.expressed_neoantigen_count / d.expressed_missense_snv_count\n",
    "\n",
    "extra_cols = ['missense_expression_ratio','neoantigen_expression_ratio', 'exonic_expression_ratio', 'missense2exonic_snv_ratio','expressed_neoantigen2missense_ratio']\n",
    "## create recentered versions of ratios\n",
    "for col in extra_cols:\n",
    "    col_centered = '{}_centered'.format(col)\n",
    "    col_rescaled = '{}_rescaled'.format(col)\n",
    "    d[col_centered] = d[col] - np.mean(d[col])\n",
    "    d[col_rescaled] = d[col_centered]/np.std(d[col_centered])\n",
    "\n",
    "## append extra-cols to key var lists\n",
    "vars_centered.extend(['{}_centered'.format(col) for col in extra_cols])\n",
    "vars_rescaled.extend(['{}_rescaled'.format(col) for col in extra_cols])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## center variables by mean within PD-L1 group"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:22.220909",
     "start_time": "2016-08-17T07:59:22.199735"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## identify list of variables to center\n",
    "metrics = list(cols)\n",
    "metrics.extend(extra_cols)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "metrics2 = list(metrics)\n",
    "metrics2.extend(['pd_l1'])\n",
    "assert(not 'pd_l1' in metrics)\n",
    "log_metrics2 = ['log_{}'.format(var) for var in metrics2]\n",
    "metrics2.extend(log_metrics2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:22.242582",
     "start_time": "2016-08-17T07:59:22.222780"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "grp_metrics = [var for var in metrics2 if var in d.columns]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:22.275217",
     "start_time": "2016-08-17T07:59:22.244289"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# center variables by group\n",
    "bygrp = d.loc[:, grp_metrics]\n",
    "bygrp = bygrp.groupby('pd_l1').transform(lambda x: x - x.mean())\n",
    "bygrp['patient_id'] = d.patient_id"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:22.301832",
     "start_time": "2016-08-17T07:59:22.276693"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# merge recentered variables back into original dataframe\n",
    "df = pd.merge(d, bygrp, on = 'patient_id', suffixes = ['', '_centered_by_pd_l1'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# prep data for survival analysis (long format)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:22.420738",
     "start_time": "2016-08-17T07:59:22.356618"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## prep dflong_pfs which will be used for survival analysis using stan\n",
    "df_long_pfs = survivalstan.prep_data_long_surv(df = df,\n",
    "                                               event_col = 'is_progressed_or_deceased',\n",
    "                                               time_col = 'pfs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# multivariate analysis using varying-coef model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T08:27:09.524216",
     "start_time": "2016-08-17T08:27:09.501719"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "../utils/stan/logistic_model.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_custom_coefs_expressed_missense_and_neoant_mutations.stan\n",
      "../utils/stan/logistic_model_by_group.stan\n",
      "../utils/stan/pem_survival_model_unstructured_varcoef.stan\n",
      "../utils/stan/pem_survival_model_unstructured_varcoef_hsprior.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_tvc.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_alt.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_bspline_est_xi.stan\n",
      "../utils/stan/pem_survival_model_varying_coefs3.stan\n",
      "../utils/stan/pem_survival_model_randomwalk.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_custom_coefs_rate_only.stan\n",
      "../utils/stan/pem_survival_model_gamma.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_custom_coefs_missense_and_neoant_rates.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_tvc2.stan\n",
      "../utils/stan/pem_survival_model_varying_coefs2.stan\n",
      "../utils/stan/pem_survival_model_varying_coefs4.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_bspline_est_xi2.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_bspline.stan\n",
      "../utils/stan/pem_survival_model_randomwalk2.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_bspline2.stan\n"
     ]
    }
   ],
   "source": [
    "models = survivalstan.utils.read_files('../utils/stan')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "survstan_pfs_varcoef = functools.partial(\n",
    "    survivalstan.fit_stan_survival_model,\n",
    "    df = df_long_pfs,\n",
    "    model_code = models['pem_survival_model_unstructured_varcoef.stan'],\n",
    "    timepoint_end_col = 'end_time',\n",
    "    event_col = 'end_failure',\n",
    "    sample_col = 'patient_id',\n",
    "    chains = 4,\n",
    "    iter = 10000,\n",
    "    grp_coef_type = 'vector-of-vectors',\n",
    "    seed = seed,\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## summarize varying-coef model for missense_snv_count by liver-mets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NOT reusing model.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ran in 91.687 sec.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/tavi/miniconda2/lib/python2.7/site-packages/stanity/psis.py:228: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  elif sort == 'in-place':\n",
      "/home/tavi/miniconda2/lib/python2.7/site-packages/stanity/psis.py:246: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n",
      "  bs /= 3 * x[sort[np.floor(n/4 + 0.5) - 1]]\n"
     ]
    }
   ],
   "source": [
    "livermet = survstan_pfs_varcoef(formula = 'log_missense_snv_count_centered', \n",
    "                               group_col = 'liver_mets',\n",
    "                              )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{{{hr_missense_by_livermet_pfs}}}\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhwAAAGHCAYAAAD7t4thAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8TXfi//H3tUQkEUtSWrGkBNeSTQhtKSOYohVqqtqM\nRqnW3h8lpYN21FJRLYpSbTqkaNTSRMswtVR1mYihquiUKI09loglQnp/f5icb67skpOLvJ6PRx7i\n3M8553M+9yR538/5nPOx2Gw2mwAAAExUxtEVAAAA9z4CBwAAMB2BAwAAmI7AAQAATEfgAAAApivn\n6AoAZktLS9PevXt13333qWzZso6uDgDcFTIyMnTmzBk1a9ZMzs7ORd4egQP3vL179yosLMzR1QCA\nu9LSpUvVokWLIm+HwIF73n333Sfp5g/N/fff7+DaAMDd4eTJkwoLCzN+hxYVgQP3vMzLKPfff79q\n1arl4NoAwN2luC5FM2gUAACYjsABAABMR+AAAACmI3AAgIN4e3vL29vb0dUASgSBAwAAmI7AAQAA\nTEfgAAAApiNwAAAA0xE4AACA6XjSKAA4yG+//eboKgAlhh4OAABgOgIHAAAwHYEDAACYjjEcAO4Y\nERERSk5OLvR6ly5dkiS5ubnd1n49PT0VGRl5W+sCKBgCB4A7RnJyss6cPi2XQq539X//Wq5cKfQ+\nC78GgNtB4ABwR3GR9FS5wv1q+uzGDek21su6riNkzqPC3SooDRjDAQAATEfgAAAApiNwAAAA0xE4\nAACA6QgcAADAdNylAgAOwt0pKE3o4QAAAKYjcAAAANMROAAAgOkIHAAAwHQEDgAAYDoCBwA4iLe3\ntzGfCnCvI3AAAADTETgAAIDpCBwAAMB0BA4AAGA6AgcAADAdc6kAgIMwlwpKE3o4AACA6QgcAADA\ndAQOAABgOgIHAAAwHYEDAACYjsABAA7CXCooTQgcAADAdAQOAABgOgIHAAAwHYEDAACYjsABAABM\nx1wqAOAgzKWC0oQeDgCFEhUVpaioKEdX455E2+JeRuAAUCjbt2/X9u3bHV2NexJti3sZgQMAAJiO\nwAEAAExH4AAAAKYjcACAgzCXCkoTAgcAADAdgQMAAJiOwAEAAExH4AAAAKYjcAAAANMxlwoAOAhz\nqaA0oYcDAACYjsABAABMR+AAAACmI3AAAADTETgAAIDpCBwA4CDMpYLShMABAABMR+AAAACmI3AA\nAADTETgAAIDpCBwAAMB0zKUCAA7CXCooTejhAAAApiNwAAAA0xE4AACA6UwNHGvWrJHValVwcLBS\nU1PtXsvIyJDVatXcuXOLdV+///57rmXGjRunkJCQYtlfcejbt6+sVqueeeaZHF8fN26crFar2rdv\nX+htHzt2THPnzlVSUlIRa5m7xYsX61//+pdp21+zZo0aN26s48ePm7YPAEDJKJEejtTUVC1atMj0\n/VgsljxfHzJkiObNm2d6PQrDzc1Nu3fvzhaU0tLStGHDBrm5ud3WdjMDR14BrKjMDhzt27dXTEyM\n7rvvPtP2AQAoGSUSOB555BFFR0fr3LlzJbG7XNWuXVtWq7VE95menp7n640aNVKdOnUUGxtrt3zD\nhg2yWCxq06bNbe3XZrPlG8DudFWrVpWfn5/Kly/v6KoApmAuFZQmpgcOi8WiwYMHS5Lmz5+fb/k9\ne/aoX79+CgwMVGBgoPr166c9e/YUS13Gjh2rDh06SLoZBFq1aqXp06dnK7du3TpZrVYdOHDAWBYf\nH69+/fqpefPmCgwM1IABA/Trr7/arde3b189++yz2rJli3r27Ck/Pz8tX74833qFhoZq7dq1dsvi\n4uLUuXNnVaxYMVv5jIwMLVy4UF26dJGvr6/atm2r6dOnG+EmPj5e4eHhkqTnn39eVqtVjRs31o4d\nO4zjCw8P10MPPaTAwED17NlTn3/+ebb9LF68WF27dpW/v7+Cg4PVq1cvffXVV5KkDh066MSJE4qL\ni5PVapXVatW4ceMkSUePHlVERIRCQkLk7++vjh076o033tDFixfttr9nzx71799frVq1MspNmjTJ\neH316tWyWq12l1TWrl2rnj17KjAwUEFBQXriiSe0YsWKfNsYAOBYJfIcjurVqyssLExLlizRgAED\n9MADD+RY7sCBA+rbt698fHwUGRkpSVq4cKH69u2rFStWqFGjRkWqh8ViMT71Ozk56bHHHtMXX3yh\niIgIu96AuLg4NWzY0OgN2bp1q4YOHao//elPevvttyVJH3zwgcLCwrR27VrVqFHDWPe3337TlClT\nNGTIENWuXVuVK1fOt17du3fXe++9p927dysgIECnTp3S999/r6ioqGw9H5I0evRobd26VS+++KIC\nAgKUmJioWbNm6dixY5ozZ46aNGmiiRMn6s0339SECRPk6+srSapfv76km4GgU6dOGjhwoMqWLauE\nhASNHz9e165d09NPP220QWRkpIYNG6agoCClpaXpl19+UUpKiqSb4fGFF15Q48aNNXz4cEk3eyQk\n6fTp06pRo4bGjRunKlWqKCkpSQsWLNCLL76oTz/9VJJ05coVDRw4UP7+/oqMjJSLi4uOHTum//zn\nPzm+X5KUkJCgiIgIhYeHKyIiQjabTYmJidmCDADgzlNiD/4aOHCgYmJiNHfuXE2ZMiXHMvPnz1eF\nChW0ePFiY+zCQw89pJCQEM2bN09z5swp1jqFhoYqJiZG3333nR555BFJ0rlz57R9+3aNGjXKKDd1\n6lS1atXKboBrq1atFBISoqioKOOTvSRduHBBH3/8caHCUa1atRQUFKTPP/9cAQEBiouL0/3336/W\nrVtnCxwJCQlav369IiMj1b17d0k328jd3V0RERE6cOCArFarfHx8ZLPZVK9ePfn5+dltY9CgQcb3\nNptNwcHBOn36tJYvX24Ejh9//FGNGjUyeqck6dFHHzW+t1qtcnJyMi57ZNWiRQu1aNHC+H9gYKBq\n166tsLAwo36ZQWH06NFq2LChJKlly5bq0aNHru20Z88eubu7a+zYscayhx9+OO/GBQDcEUrsttjK\nlSvr+eefV2xsbK5P10tISFD79u3tBkq6ubmpQ4cOio+PL/Y6NW/ePNv4iS+//FI2m02PP/64JOnI\nkSM6evSoHn/8cWVkZBhfFSpUUEBAgHGZIpOXl9dt9cSEhoZq/fr1Sk9PV1xcnJ544okcy33zzTdy\ncnLSn//8Z7v6PPLII7LZbEpISMh3X0eOHNGoUaP06KOPqmnTpmratKk+++wzHT582Cjj6+urAwcO\naPLkyfr++++VlpZW4GO5fv26FixYoC5dusjf319NmzZVWFiYLBaLEhMTJd28du3u7q6JEycqLi5O\nJ0+ezHe7vr6+unjxosaMGaOtW7dmu/MJAHDnKtFHm/fr10+ffPKJ5syZoxkzZmR7PSUlJcc7Ejw9\nPU3rNu/evbuioqKUlpYmZ2dnxcXFqXXr1qpevbok6ezZs5Kkv/3tb3rttdfs1rVYLNkuD93uHRVd\nunTRlClTNG/ePB08eDDX3pxz584pPT1d/v7+2V6zWCy6cOFCnvu5cuWKnn/+ebm4uGjMmDGqXbu2\nypcvr2XLlmn16tVGuR49eig9PV0rV67U8uXLVbZsWbVr105jx46Vl5dXnvuYOXOmli5dqmHDhikg\nIECurq46efKkhg0bZowzcXNz0+LFizV//nxNmjRJly5dUoMGDTR8+HB17tw5x+22bNlSs2fPVnR0\ntIYNG2YsGzt2bJEvt6HgLl26pLS0NPXv37/Yt52cnKyyxb7VvKVLupqcbMrx5Cezdy9z38nJyXJ2\ndi7xegAloUQDh4uLi1588UVFRkbm+MNduXJlJScnZ1uenJwsd3d3U+oUGhqquXPnauPGjfLz89NP\nP/1kjB+RpCpVqkiSRo0alWP3/a13UNzunSGZPTmLFi2Sr6+vHnzwwRzLValSRc7Ozlq2bJlsNlu2\n1zODUm512b17t06cOKFly5YpMDDQWH7jxo1sZXv37q3evXsrNTVV27dv11tvvaVRo0YpJiYmz2NZ\nt26devbsqZdeeslYdvny5WzlrFar5syZoz/++EN79+7VwoULNXLkSMXGxsrHxyfHbXfu3FmdO3fW\n1atXFR8frxkzZmjgwIHatm1bnnUC7kS1atVydBWAElPik7c9++yzWrx4sWbNmpXtD2LLli319ddf\n68qVK3JxcZF089PU5s2b1bp1a1PqU7t2bQUGBio2NlaHDx+Wi4uLOnXqZLxer149eXl56eDBgxo4\ncKApdcgUFham9PT0XC+nSFLbtm314Ycf6uLFi3m2iZOTk2w2m65du2a3/OrVq5KksmX/73NkSkqK\nNm/enOu2KlWqpC5duujHH3+0CxtOTk45XmpJS0uz274krVq1KtcwVqZMGfn5+WnEiBHatGmTDh06\nlGvgyFSxYkW1a9dOR48e1dSpU3X+/Hlj0CrM5ebmJjc3N0VFRRX7tvv376/Lp08X+3bz4iTJ1dPT\nlOMpLEf0sgAlpcQDh5OTk4YMGaIJEyZk+wM0ZMgQff311woPDzf+uC9atEjXrl3T0KFD8922zWbT\ntm3b5Onpabe8UqVKeQ4uDA0N1aRJk/TLL7+oU6dO2W5FnThxooYOHar09HR16dJFVatWVXJysnbt\n2qWaNWuqX79+BTz6vAUFBSkoKCjPMsHBweratatefvllhYeHy8/PT2XKlFFSUpK2bdumMWPGqG7d\nuvL29la5cuW0atUqubu7y8nJSfXq1VNgYKBcXV01adIkDR8+XJcvX9aCBQtUrVo1Xbp0ye6YXV1d\nFRAQIA8PDx0+fFixsbFq27atUcbHx0c7d+7U1q1b5enpqapVq8rLy0tt27bV559/rgYNGqhu3bra\nuHGjdu/ebXccW7duVUxMjDp27KhatWrpypUrio6OlpubmwICAnI89jlz5ig5Odm45HXixAlFR0er\ncePGhA0AuMM5ZHr6J598Uh9++KGOHj1qt7xRo0ZasmSJZs2apbFjx8pmsykwMFCffPKJca0zLxaL\nRZMnT8623MfHx3jORU6fsrt27aopU6bo3LlzCg0NzfZ6u3bttHTpUr3//vuaMGGC0tLS5OnpqYCA\nAHXr1i1bHQqjIOVvLTNz5kxFR0dr1apVWrhwoZycnOTl5aU2bdrIw8ND0s1LLxMnTtSiRYv03HPP\nKSMjQ0uWLFHLli01b948TZ8+XS+//LKqV6+u5557ThcuXLB7Cmvz5s21evVqxcXFKTU1VdWrV1eP\nHj2MsRPSzctMEydO1MiRI5WWlqYePXpo2rRpGj9+vCRp9uzZRvu98847euqpp4x169atq4oVK+r9\n99/XmTNn5OrqKl9fX0VFRdndZpyVv7+/oqOjNW3aNKWkpMjDw0Nt2rTRiBEjCtjaAABHsdhyGggA\n3EOSkpIUEhKiTZs2cc28GGR2+5t5SeWpcoX7LPTZ/8YgFXa9zHVdq1e/oy6p3Al1AYr7dyezxQIA\nANMROADAQZhLBaUJgQMAAJiOwAEAAExH4AAAAKYjcAAAANMROAAAgOkc8uAvAIBynTkbuBfRwwEA\nAExH4AAAAKYjcAAAANMROAAAgOkIHAAAwHQEDgBwEOZSQWlC4AAAAKYjcAAAANMROAAAgOkIHAAA\nwHQEDgAAYDrmUgEAB2EuFZQm9HAAAADTETgAAIDpCBwAAMB0BA4AAGA6AgcAADAdgQMAHIS5VFCa\nEDgAAIDpCBwAAMB0PPgLQKG0adPG0VW4Z9G2uJcROAAUSv/+/R1dhXsWbYt7GZdUAACA6ejhAAAH\nYS4VlCb0cAAAANMROAAAgOkIHAAAwHQEDgAAYDoCBwAAMB2BAwAchLlUUJoQOAAAgOkIHAAAwHQE\nDgAAYDoCBwAAMB2BAwAAmI65VADAQZhLBaUJPRwAAMB0BA4AAGA6AgcAADAdgQMAAJiOwAEAAExH\n4AAAB2EuFZQmBA4AAGA6AgcAADAdgQMAAJiOwAEAAExH4AAAAKZjLhUAcBDmUkFpQg8HAAAwHYED\nAACYjsABAABMR+AAAACmI3AAAADTcZcKgDvKFUmf3bhR6HV0G+tlruta6LWKR+Y8KtytgtKAwAHg\njuHp6Xlb69kuXZIkubq5FXpd1yLsF0DBETgA3DEiIyMdXQUAJmEMBwAAMB2BAwAAmI7AAQAATMcY\nDgBwEO5OQWlCDwcAADAdgQMAAJiOwAEAAExH4AAAAKYjcAAAANMROADAQby9vY35VIB7HYEDAACY\njsABAABMR+AAAACmI3AAAADTETgAAIDpmEsFAByEuVRQmtDDAQAATEfgAAAApiNwAAAA0xE4AACA\n6QgcAADAdAQOAHAQ5lJBaULgAAAApiNwAAAA0xE4AACA6QgcAADAdAQOAABgOuZSAQAHYS4VlCb0\ncAAAANMROAAAgOm4pAIA/xMREaHk5OQib+fSpUuSJDc3tyJvKytPT09FRkYW6zaBkkLgAID/SU5O\n1unTZ1S2QsUibSfj2lVJ0rWM4qiV/TaBuxWBAwCyKFuhoqoHdS/SNk7vjJOkIm8np20CdyvGcACA\ng3w6Y6Q+nTHS0dUASgSBAwAAmI7AAQAATEfgAAAApiNwAAAA0xE4AACA6bgtFgAcpM+Ydx1dBaDE\n0MMBAABMR+AAAACmI3AAAADTETgAAIDpCBwAAMB0BA4AcBDmUkFpQuAAAACmI3AAAADTETgAAIDp\nCBwAAMB0BA4AAGA65lIBAAdhLhWUJvRwAAAA0xE4AACA6QgcAADAdAQOAABgOgIHAAAwHYEDAByE\nuVRQmhA4AACA6QgcAADAdAQOAABgOgIHAAAwHYEDAACYjrlUAMBBmEsFpQk9HADuWFFRUYqKinJ0\nNWAC3tvSh8AB4I61fft2bd++3dHVgAl4b0sfAgcAADAdgQMAAJiOwAEAAExH4AAAB2EuFZQmBA4A\nAGA6AgcAADAdgQMAAJiOwAEAAExH4AAAAKZjLhUAcBDmUkFpQg8HAAAwHYEDAACYjsABAABMR+AA\nAACmI3AAAADTETgAwEGYSwWlCYEDAACYjsABAABMR+AAAACmI3AAAADTETgAAIDpmEsFAByEuVRQ\nmtDDAQAATEfgAAAApiNwAAAA0+UZOL788ktZrVYlJCTYLT979qysVqvatGmTbZ2lS5fKarXq4MGD\nkqRx48YpJCTEeP3YsWOaO3eukpKSsq3boUMHRURE3NaBrFmzRlarVb///nuuZW6ti6P17dtXVqtV\nzzzzTI6vjxs3TlarVe3bty/0tvNq5+KyePFi/etf/zJt+2vWrFHjxo11/Phx0/YBACgZeQaOFi1a\nSFK2wLFjxw5VrFhRZ8+e1eHDh+1eS0hIUNWqVeXj4yNJGjJkiObNm2e8nvmHMK9gcLssFkuer99a\nlzuBm5ubdu/ena090tLStGHDBrm5ud3Wds1s50xmB4727dsrJiZG9913n2n7AACUjDwDR40aNVSn\nTh3t2LHDbnlCQoJat26tOnXqZAsjCQkJCgoKMv5fu3ZtWa1W4/82my3fYGCWW+tSEtLT0/N8vVGj\nRqpTp45iY2Ptlm/YsEEWiyXHXqSCcGQ7F5eqVavKz89P5cuXd3RVAFMwlwpKk3zHcLRo0UK7d+/W\nH3/8YSzbsWOHWrRooebNm9uFkSNHjujMmTNq2bKlsWzs2LHq0KGDJCk+Pl7h4eGSpOeff15Wq1WN\nGzfOFmjWrVunrl27KjAwUL169dLOnTuLdpQ51CU9PV2tWrXS9OnTs5Vbt26drFarDhw4YCyLj49X\nv3791Lx5cwUGBmrAgAH69ddf7dbr27evnn32WW3ZskU9e/aUn5+fli9fnm+9QkNDtXbtWrtlcXFx\n6ty5sypWrJitfEZGhhYuXKguXbrI19dXbdu21fTp041wk187r1u3TuHh4XrooYcUGBionj176vPP\nP8+2n8WLF6tr167y9/dXcHCwevXqpa+++krSzctfJ06cUFxcnKxWq6xWq8aNGydJOnr0qCIiIhQS\nEiJ/f3917NhRb7zxhi5evGi3/T179qh///5q1aqVUW7SpEnG66tXr5bVarW7pLJ27Vr17NlTgYGB\nCgoK0hNPPKEVK1bk28YAAMfK9zkcLVu21Jo1a/Tzzz/L19dXqamp+vXXX9WiRQtVrlxZ8+fPN8rG\nx8fLYrHYBQ6LxWJ80m7SpIkmTpyoN998UxMmTJCvr68kqX79+kb5hIQEHT58WCNHjpSTk5NmzZql\nwYMHa/Pmzbd9eSGnujg5Oemxxx7TF198oYiICLvegLi4ODVs2NDoDdm6dauGDh2qP/3pT3r77bcl\nSR988IHCwsK0du1a1ahRw1j3t99+05QpUzRkyBDVrl1blStXzrde3bt313vvvafdu3crICBAp06d\n0vfff6+oqKhsPR+SNHr0aG3dulUvvviiAgIClJiYqFmzZunYsWOaM2dOvu189OhRderUSQMHDlTZ\nsmWVkJCg8ePH69q1a3r66aeNNoiMjNSwYcMUFBSktLQ0/fLLL0pJSZEkzZ8/Xy+88IIaN26s4cOH\nS7rZIyFJp0+fVo0aNTRu3DhVqVJFSUlJWrBggV588UV9+umnkqQrV65o4MCB8vf3V2RkpFxcXHTs\n2DH95z//yfH9km6eGxEREQoPD1dERIRsNpsSExOzBRkAwJ2nQIHDZrMpISFBvr6+2rFjhypUqKCm\nTZuqcuXKOn78uI4fP66aNWsqISFBbm5uaty4cY7bcnNzk4+Pj2w2m+rVqyc/P79sZS5fvqy4uDgj\nXHh4eOgvf/mLvv76a3Xr1q2Ih2svNDRUMTEx+u677/TII49Iks6dO6ft27dr1KhRRrmpU6eqVatW\nmjt3rrGsVatWCgkJUVRUlPHJXpIuXLigjz/+WI0aNSpwPWrVqqWgoCB9/vnnCggIUFxcnO6//361\nbt06W+BISEjQ+vXrFRkZqe7du0uSHnroIbm7uysiIkIHDhyQ1WrNs50HDRpkfG+z2RQcHKzTp09r\n+fLlRuD48ccf1ahRIw0ePNgo++ijjxrfW61WOTk5GZc9smrRooUx/keSAgMDVbt2bYWFhRn1ywwK\no0ePVsOGDSXdPNd69OiRazvt2bNH7u7uGjt2rLHs4YcfzrtxAQB3hHwDR61atXT//fdrx44dev75\n55WQkCA/Pz+VK1dO3t7e8vDw0I4dOxQaGqqEhAQ1b968SGMHAgIC7HoyMv8YmXGnQvPmzY3xE5mB\n48svv5TNZtPjjz8u6eZloqNHj2rQoEHKyMgw1q1QoYICAgKyXQ7y8vIqVNjIFBoaqrfffluvvfaa\n4uLi9MQTT+RY7ptvvpGTk5P+/Oc/29XnkUceMYJhfuNUjhw5otmzZyshIUHJycnG5bIKFSoYZXx9\nfbV8+XJNnjxZISEhCgwMlLOzc4GO5fr16/roo48UGxur48eP69q1a5Ju9lgkJibKarXK29tb7u7u\nmjhxop599lkFBwfr/vvvz3O7vr6+unjxosaMGaNu3bopKChIlSpVKlCdcHe6dOmS0tLS1L9//xLZ\nX3JysmyWO/NpAX/cSFdycnKJtYXZkpOTC/w7BfeGAj3avEWLFvrmm28k3fyE3bZtW+O1oKAgYxDp\nsWPH1KdPnyJV6NZLEE5OTpJk/NEqbt27d1dUVJTS0tLk7OysuLg4tW7dWtWrV5d08xZgSfrb3/6m\n1157zW5di8WiBx54wG7Z7d5R0aVLF02ZMkXz5s3TwYMHNWfOnBzLnTt3Tunp6fL398/2msVi0YUL\nF/Lcz5UrV/T888/LxcVFY8aMUe3atVW+fHktW7ZMq1evNsr16NFD6enpWrlypZYvX66yZcuqXbt2\nGjt2rLy8vPLcx8yZM7V06VINGzZMAQEBcnV11cmTJzVs2DBjnImbm5sWL16s+fPna9KkSbp06ZIa\nNGig4cOHq3Pnzjlut2XLlpo9e7aio6M1bNgwY9nYsWNvK+QBAEpOgQJHcHCwvvzyS+3evVv79u3T\nyJH/N6o6KChIy5cvz3H8Rm7upLsnQkNDNXfuXG3cuFF+fn766aefFBkZabxepUoVSdKoUaNy7L6/\n9Q6K2z02Nzc3dejQQYsWLZKvr68efPDBHMtVqVJFzs7OWrZsmWw2W7bXM4NSbnXZvXu3Tpw4oWXL\nlikwMNBYfuPGjWxle/furd69eys1NVXbt2/XW2+9pVGjRikmJibPY1m3bp169uypl156yVh2+fLl\nbOWsVqvmzJmjP/74Q3v37tXChQs1cuRIxcbGGrdV36pz587q3Lmzrl69qvj4eM2YMUMDBw7Utm3b\n8qwT7k5ubm5yc3NTVFRUieyvf//+OptyqUT2JRVuLpUy5ZzkUbnk2sJs90pPDQquwD0cNptNH3zw\ngaSblz0yBQUFadq0aVq/fr2cnZ2NAYq5cXJyks1mM63HorBq166twMBAxcbG6vDhw3JxcVGnTp2M\n1+vVqycvLy8dPHhQAwcONLUuYWFhSk9Pz/VyiiS1bdtWH374oS5evKjWrVvnWi63dr569aokqWzZ\nssaylJQUbd68OddtVapUSV26dNGPP/5oFzacnJyUlpaWrXxaWprd9iVp1apVuYaxMmXKyM/PTyNG\njNCmTZt06NChXANHpooVK6pdu3Y6evSopk6dqvPnzxuDVgEAd54CBY569erJw8NDW7ZsUbNmzexu\n1WzSpIlcXFy0ZcsWtW7dOtsfmlt5e3urXLlyWrVqldzd3eXk5KR69erJxcWlaEeimwMgt23bJk9P\nT7vllSpVynNwYWhoqCZNmqRffvlFnTp1ynYr6sSJEzV06FClp6erS5cuqlq1qpKTk7Vr1y7VrFlT\n/fr1K3LdpZvhLeszTHISHBysrl276uWXX1Z4eLj8/PxUpkwZJSUladu2bRozZozq1q2bazsHBgbK\n1dVVkyZN0vDhw3X58mUtWLBA1apV06VL//fJbuLEiXJ1dVVAQIA8PDx0+PBhxcbG2l1O8/Hx0c6d\nO7V161Z5enqqatWq8vLyUtu2bfX555+rQYMGqlu3rjZu3Kjdu3fbHcfWrVsVExOjjh07qlatWrpy\n5Yqio6NLmI7qAAAaeklEQVTl5uZmF2izmjNnjpKTk41LXidOnFB0dLQaN25M2ACAO1yBp6dv0aKF\nNm7caHf3gXTz02lgYKC+++67XC+nZP1kW6VKFU2cOFGLFi3Sc889p4yMDC1ZskQtW7bMdhtk1vUL\ncqnCYrFo8uTJ2Zb7+PgYz7nIaTtdu3bVlClTdO7cOYWGhmZ7vV27dlq6dKnef/99TZgwQWlpafL0\n9FRAQEC2O2cKe0mloMeV1cyZMxUdHa1Vq1Zp4cKFcnJykpeXl9q0aSMPDw9JebfzvHnzNH36dL38\n8suqXr26nnvuOV24cMHuKazNmzfX6tWrFRcXp9TUVFWvXl09evQwxk5INy8zTZw4USNHjlRaWpp6\n9OihadOmafz48ZKk2bNnG+33zjvv6KmnnjLWrVu3ripWrKj3339fZ86ckaurq3x9fRUVFWV3m3FW\n/v7+io6O1rRp05SSkiIPDw+1adNGI0aMKGBrAwAcxWLLaSAAcA9JSkpSSEiINm3apFq1ajm6OiiE\nzOv8JT2Go3pQ9yJt5/TOOEkq8nZu3ea9OIbjXjmee1Fx/+68M+//AgAA9xQCBwA4CHOpoDQhcAAA\nANMROAAAgOkIHAAAwHQEDgAAYDoCBwAAMF2BH/wFAChehZlLBbjb0cMBAABMR+AAAACmI3AAAADT\nETgAAIDpCBwAAMB0BA4AcBDmUkFpQuAAAACmI3AAAADTETgAAIDpCBwAAMB0BA4AAGA65lIBAAdh\nLhWUJvRwAAAA0xE4AACA6QgcAADAdAQOAABgOgIHAAAwHYEDAByEuVRQmhA4AACA6QgcAADAdDz4\nC8Adq02bNo6uAkzCe1v6EDgA3LH69+/v6CrAJLy3pQ+XVAAAgOno4QAAB2EuFZQm9HAAAADTETgA\nAIDpCBwAAMB0BA4AAGA6AgcAADAdgQMAHIS5VFCaEDgAAIDpCBwAAMB0BA4AAGA6AgcAADAdgQMA\nAJiOuVQAwEGYSwWlCT0cAADAdAQOAABgOgIHAAAwHYEDAACYjsABAABMR+AAAAdhLhWUJgQOAABg\nOgIHAAAwHYEDAACYjsABAABMR+AAAACmYy4VAHAQ5lJBaUIPBwAAMB2BAwAAmI7AAQAATEfgAAAA\npmPQKABkkXHtqk7vjCvyNiQVeTvZt+lWbNsDShqBAwD+x9PTs1i2c+nSzX/d3PIOCJ999pkk6amn\nnirAVt2KrX6AIxA4AOB/IiMjS3R/mzdvliRFRUWV6H4BR2AMBwAAMB2BAwAAmI7AAQAATEfgAAAA\npmPQKAA4yG+//eboKgAlhh4OAABgOgIHAAAwHYEDAACYjsABAABMR+AAAACmI3AAgIN4e3vL29vb\n0dUASgSBAwAAmI7AAQAATEfgAAAApuNJo7jnZWRkSJJOnjzp4JoAOUtKSnJ0FYBsMn9nZv4OLSoC\nB+55Z86ckSSFhYU5uCaAvQoVKkiSQkJCHFwTIHdnzpxR3bp1i7wdi81msxVDfYA7Vlpamvbu3av7\n7rtPZcuWdXR1AOCukJGRoTNnzqhZs2ZydnYu8vYIHAAAwHQMGgUAAKYjcAAAANMROAAAgOkIHAAA\nwHTcFou71smTJzV16lR99913stlsevjhh/Xaa6/pgQceyHfd9PR0vfvuu1q7dq1SU1PVuHFjjR49\nWi1atCiBmt+ZitKeVqs12zKLxaI1a9bk+FppcOrUKX3wwQf6+eefdeDAAaWlpWnz5s2qWbNmvuty\nftorSltybtr75z//qbVr1+rnn3/W+fPn9cADD6hz58566aWX5Orqmue6RT0vuUsFd6W0tDR1795d\nFSpU0MiRIyVJ7777rq5du6a4uLh8b+F65ZVX9M033ygiIkK1atXS0qVLtW3bNsXExJTKX0JFbU+r\n1apevXrp6aeftlveqFEj41kTpU18fLxGjRqlpk2bKiMjQ99++602bdpUoD+SnJ/2itKWnJv2nn76\nadWoUUOdOnXS/fffr/379+u9995T/fr19emnn+a5bpHPSxtwF/rHP/5ha9Kkie3o0aPGst9//93W\npEkT28cff5znuvv377c1atTItmbNGmPZjRs3bH/+859tgwcPNqvKd7SitKfNZrM1atTINmvWLBNr\neHdbsWKFzWq12o4dO5ZvWc7PvBWmLW02zs1bnTt3LtuyNWvW2KxWq+2HH37Idb3iOC8Zw4G70pYt\nW+Tv76/atWsby2rVqqXmzZtr06ZNea67adMmlS9fXl26dDGWlS1bVt26ddP27dt1/fp10+p9pypK\ne6J4cX7CTFWrVs22zNfXVzabTadOncp1veI4LwkcuCsdPHhQDRo0yLbcx8dHhw4dynPdQ4cOqVat\nWtm6U318fHT9+nUdPXq0WOt6NyhKe2Zavny5fH19FRAQoPDwcCUkJBR3NUsFzs/ix7mZt/j4eFks\nFtWvXz/XMsVxXjJoFHelCxcuqHLlytmWV65cWRcvXsxz3ZSUlBzXrVKlirHt0qYo7SlJoaGhat++\nvapXr67jx4/ro48+Ur9+/fTxxx+rZcuWZlT5nsX5Wbw4N/N26tQpvffee3r44YfVtGnTXMsVx3lJ\n4ABQZNOnTze+DwoKUocOHfTEE09o9uzZ+uSTTxxYM5R2nJu5u3LligYPHqzy5ctr6tSppu+PSyq4\nK1WuXFkpKSnZlqekpMjd3T3Pdd3d3XNcNzOhZyb20qQo7ZkTV1dXtWvXTj/99FNxVK9U4fw0F+fm\nTdeuXdNLL72kY8eO6aOPPlKNGjXyLF8c5yWBA3clHx8fHTx4MNvygwcP5nkdMnPdpKQkXbt2Ldu6\n5cuXV506dYq1rneDorQnihfnJ8x248YNDR8+XPv27dOiRYvk4+OT7zrFcV4SOHBX6tChg3788Ucl\nJSUZy5KSkrRr1y6FhITku+7169e1fv16Y1lGRobWr1+vNm3aqHz58qbV+05VlPbMyaVLl7R161b5\n+fkVZzVLBc5Pc5X2c9Nms+mVV15RfHy85s+fX+B2KI7zsuwbb7zxxu1WHHCURo0aad26ddqwYYOq\nV6+uw4cP6/XXX1fFihU1efJk4+Q/fvy4WrVqJYvFYgwQu++++5SYmKhly5apSpUqunjxot5++239\n9NNPevvtt+Xp6enIQ3OIorRnVFSUPv/8c129elXnz59XfHy8xo8fr+PHj2v69OkFelLpvWrDhg06\ndOiQdu7cqb1798rb21vHjx/X+fPn5eXlxflZCLfTlpyb2b3xxhuKjY3VwIED5ePjo1OnThlfFotF\nbm5upp2XDBrFXalixYpavHixpk6dqldffdV4FPe4ceNUsWJFo5zNZjO+snrrrbf07rvvavbs2UpN\nTZXVatVHH31UKp/iKBWtPR988EF99dVX2rhxo1JTU+Xm5qagoCBNmzZNzZo1c8Th3DFefvllWSwW\nSTcfpz1p0iRJUsuWLbVkyRLOz0K4nbbk3Mzum2++kcVi0YIFC7RgwQK714YOHaphw4aZdl7yaHMA\nAGA6xnAAAADTETgAAIDpCBwAAMB0BA4AAGA6AgcAADAdgQMAAJiOwAEAAEzHg7+AYjR37lzNnTtX\nkjRs2DANGzbM7vXMB+RYLBbt378/2/JMFotFrq6uatCggXr37q2ePXsWaP+//fab3nrrLe3Zs0fn\nz5+XzWbTa6+9pueee64oh1UgWY89U5kyZVS5cmX5+/vrhRdeUIsWLUyvx93qvffe07x58yRJ0dHR\n9+zU6ceOHbN7XH7VqlW1bds2u0dj7927V3/5y1+M/3t5eWnTpk2m1CfznPXy8irwz1lO+vbtqx07\ndtj9bMfHxxs/e1l/H8THxys+Pl6S9OSTT6pmzZpFOYS7BoEDMEHmExEL89qtyy9fvqxdu3Zp165d\nunz5sv7617/mu9+IiAjt2bPH2FaZMiXfiZn1OGw2my5cuKCtW7dq27ZtmjNnjjp27FjidbobWCwW\n46s0yDzOCxcu6J///KeeeOIJ47Vly5bZlTFTZuAIDg4uUuDITU7va3x8vObOnSuLxaJWrVqVmsDB\nJRWghOX2cN/M5fv379fu3bs1dOhQSTd/YUVHRxdo2/v27ZPFYlG9evX0448/at++fcXau3H9+vUC\nlRs6dKj279+vhIQEPf3005JuHt/06dPzXTc9Pb1IdbxdGRkZ+uOPPxyyb+nmJ+D9+/dr3759d33v\nRmHeQ5vNpuXLlxv/T01N1fr1640/0GY9DDvruWxWsAkODjbe08yf59KMwAHcgSpUqKB+/fpJuvkL\n98SJE3mWX7NmjaxWqzIyMmSz2XTo0CH5+fnJarVqx44dkqTz589r6tSp6ty5s3x9fdW8eXP16dNH\nq1evtttWfHy8rFarrFar5syZowULFqhDhw5q0qSJdu/eXajjcHV11ciRI43jSEpK0oULFyTdnH3S\narUqJCRECQkJ6tOnj/z9/fX6668b669evVrPPPOMmjdvLl9fX3Xq1ElTp07V+fPn7fZz48YNzZgx\nQ23atFFgYKBeeOEFHTlyxDiOrF34mW1ltVr16aef6q233lKbNm3UrFkznTx5UpJ06tQpvf766woJ\nCVGzZs0UHBysgQMHKiEhwW6/58+f19///nd17NhRAQEBCgoK0mOPPaZXXnlFv/32m1Fu48aNCgsL\n00MPPSRfX1+1adNGf/3rX/Xxxx8bZebOnWvUK/M9k24GoX/84x968sknFRgYKD8/P3Xr1k1z5szR\n1atX7eqTuX7fvn319ddfq1evXvL391enTp304YcfFvh9K2i7F+Q9zEuNGjVUrlw57dq1S7/++quk\nm+/P1atX5eXllWvYKI5z+YsvvpDVapXFYpHNZrMrmxnSv//+e7300kvq0KGDAgMD1axZM7Vv315j\nxozR0aNH8z2+rNvM7Enp0KGD0bths9nUt29fo8z27dvVpk0bWa1WPf7443bb+vXXX41yBW3fOw2X\nVIC7gIeHR75lsn5KyzrJlSQlJyerd+/eOn78uLHsxo0b2r17t/GVORlW1m0sW7ZMKSkp2bZfGHn1\nGlgsFp07d04DBgwwPhVn7mfixIlasWKF3X6TkpK0ZMkSbdq0SStWrDDa5fXXX9eqVauMst9++636\n9u2b76WtWbNmZTu+xMREPfvss7pw4YKxLDU1Vd98842+/fZbzZw5U126dJEkvfrqq9q2bZvdfo4c\nOaIjR46oe/fu8vb21p49e/T//t//s/vjefbsWZ09e1ZpaWl6/vnns9Ura9sNGjTImHArU2JioubP\nn6+vv/5aS5culbOzs936+/fv16BBg4xlv//+u2bOnKkaNWrYXbrISWHaPXN/t76HBeXp6Sk/Pz9t\n3LhRy5cv18SJE/Xpp5/KYrHo6aef1syZM7OtU1zn8q2XOXL6/qefftK2bdvstnXq1CmtXbtW33//\nvb744gtVqVIl3+O89TzM7We1QoUK6tOnj+bOnatDhw4pISHBGPf0xRdfGOV69+6d7z7vRPRwACbJ\n+ok186ugf7TT0tKMT78WiyXbp51b9ezZU/v375fNZjOmlM7aPT9r1izjF/STTz6pf//734qNjTWu\nHX/22Wc59l6kpKRo/PjxSkhI0JYtW9SwYcNCtcGlS5c0a9Ys4/+1a9fO9gs6LS1NwcHB2rRpk3bt\n2qVBgwbpP//5j/FHr2bNmoqNjVV8fLyefPJJSdLx48c1e/ZsSTcHymaGDXd3d8XExOjf//63AgMD\nc5zxMpPNZtPVq1f1zjvvaNeuXdqwYYOqVaumKVOm6MKFC3J3d9eSJUu0Z88ebdiwQfXq1ZPNZtOb\nb76pGzduSJISEhJksVjUqVMnJSQkaOfOnYqLi9Orr76qGjVqSJJ27txphK6YmBjt3btXX3/9tRYs\nWJDv+/rFF18YYaNx48b66quv9O233+qRRx6RdPMS2pIlS7Ktd/nyZQ0aNEg7duzQ+PHjjeWxsbF5\n7q8w7Z7Vre/h4MGD89xPVn369JEkxcXFaevWrUpMTJSzs7MxnuLWn5niOpf/9Kc/5fgzs3//fi1e\nvFiS1KZNG33yySf69ttv9fPPP+vf//63XnrpJUk3Q2NcXFyBjzPT5s2bNXToUGO/0dHRdj+rzzzz\njJycnCTJ7lJT5mUmq9Wqpk2bFnq/dwICB2CSrJ+iCjIYMLOL1Wq1KiAgQPPmzVO5cuXUu3dvvfzy\ny0Wqy9dff218/+qrr8rd3V0NGzY0LtvcWibTww8/rLCwMLm6uqpGjRqqXLlygfaXGbZatGihmJgY\nSTcHsEZERNiVywwD06ZNU82aNeXs7Kw6derY1eW5555Tw4YNValSJY0dO9Zox8xPnj/88INRtkeP\nHvLz85O7u7tGjRolKe9BuqGhoerSpYucnZ1Vu3ZtWSwW/fDDD7JYLLp48aL69u0rX19fde7cWYmJ\nibLZbDp//rz27dsnSapVq5ZsNpt27dql+fPna8OGDUpPT1d4eLhx51GtWrWMfS5cuFCLFy/Wvn37\n5Ovra9f+OcnaDkOGDJGXl5eqVaum0aNH51gmk4eHh0aMGCE3Nze7gZDHjx8v8P7ya/dMub2HBfXw\nww+rbt26unz5sl599VVZLBZ169ZN7u7u+dbR7HO5evXqWrt2rfr06aOAgAAFBwfbTel++PDhAh9n\nbm4NxB4eHurWrZtsNpv+9a9/6fz58/rpp5+MSzhPPfVUkffpKFxSAUwydOjQXG+Lzcutd3lcuXKl\nyHXJvPbu4uJi94s86+j4s2fPZluvSZMmt7W/rN3E7u7uCggI0IABA3IcDOnh4SFPT0+7ZefOncux\njpUqVZKbm5tSU1ON+mYdV/DAAw/k+H1ubj2+CxcuKCMjI8+AaLFYjH1OnjxZY8eO1eHDhxUVFWX8\n8ahZs6bmz58vq9WqTp06KSwsTCtXrtTmzZu1efNm2Ww2lS1bVn369NGECRNyrV/WY8vaDl5eXsb3\nOb1vderUMerv4uJiLL927Vqu+5IK1+5Z5fQeFkbv3r01Y8YMpaSkyGKx6Jlnnsm1bEmdyzabTeHh\n4Tp06FC2S5SZ73NaWlqhtllQ4eHhWrNmja5fv66VK1ca70uFChXUvXt3U/ZZEujhAO4QmV2s+/fv\n19atWxUcHKyMjAx98cUXioyMLNK2q1WrJkm6cuWKUlNTjeVZB6PmNE6kQoUKt7W/zLtU9u3bpx9+\n+EELFizI9c6LnPaRWV/J/lN5amqqLl26JIvFYtS3atWqxuunTp0yvs8cAJqXrGMfJKlKlSoqW7as\nJKlu3bpGF3vWr3379qldu3aSJD8/P61bt06bNm3SokWLNHr0aLm4uOjEiRN6++23je1OmDBBO3bs\n0IoVKzRjxgy1a9dOGRkZWrZsmX788cdc65dbO2T9Pqf3rVy52/ssWZh2z+p2z5NMTz75pJycnGSx\nWNS0adM8LxmU1Ln8yy+/GGHDx8dHW7Zs0f79+zV//vxCbed2WK1WBQcHS5JWrFhhXE557LHH5Obm\nZvr+zULgAO5ANWrUUGRkpCpWrCjp5nMJEhMTb3t77du3N76fPn26Ll68qP/+97/6xz/+kWMZR8us\ni81mU3R0tP773/8qNTVVb731lvHpMrNM69atjbJxcXH6+eeflZKSYgw4LMxtlRUqVFDr1q1ls9l0\n5MgRzZgxQ+fOndP169eVmJiojz/+WOHh4Ub5d999V1u2bFGZMmXUqlUrPfbYY6pcubLdnUU7duzQ\nokWLlJiYKG9vb3Xu3Fn+/v7GNvK6zJH1PVmwYIGSkpKUnJxsF2aK830rTLsXp6pVq2ro0KEKCQnR\nkCFDClRHqXjO5SpVqshms+n48eO6ePGisTwzeEqSk5OTnJ2ddezYMS1cuLDA285N1pD8yy+/5HiO\nhoeHG3d2ZYbnrA9DuxtxSQW4Q9WoUUPPPfecFi5cqIyMDL377rt677338l0vp19eI0aM0Lfffqvj\nx49r5cqVWrlypfGaxWIxbme8UwQGBurpp5/WihUrdOzYMbtuZIvFIi8vLw0fPlyS5O3trb/85S9a\ntWqVzp49q169ekm6ef096zoF9dprryksLEwpKSn66KOP9NFHH9m9nvVyxvr163P8A2SxWNS2bVtJ\nNz95z5w5M8c7LlxcXBQUFJRrXbp27aq1a9dq27Zt2rt3r91D0zJ7A/r27Wu3Tn7PeclLYdq9qG6t\nT+ZgzPzKFfe5HBAQoK1btyopKcnoVRg2bJgGDx6s+vXrKzExUT///LMRbL29vXOsV2Fkrd/kyZM1\nefJkSdKBAweM5R06dFCdOnX0+++/S7rZ43a3P6mXHg6gmOV3K2ZO4wNyWz5w4EDjro6vvvpKe/bs\nyXffOW3H09NTq1atUnh4uOrWrSsnJye5uroqICBA06ZNy3Zf/+3eAlvY9fIaK/H3v/9d06ZNU0BA\ngFxdXVW+fHnVqVNH4eHhWrlypV23+RtvvKEBAwbIw8NDFStW1KOPPqo5c+YY+7j1zpi89lu/fn3F\nxsbqmWeeUZ06deTk5CR3d3c1aNBATz31lP7+978bZf/617/qoYceUo0aNYxPwQ0aNNCIESM0ZswY\nSVLTpk3Vq1cv+fj4yN3dXeXKlVO1atXUoUMHLVmyJFswylqvMmXK6P3339err76qJk2aqGLFiqpQ\noYJ8fHw0dOhQffLJJ9luiS3M+VXUds+vLXNTmIHUt5Yr7nN5/Pjxat++vSpXrmy3v7Jly2rBggV6\n9NFH5ebmpmrVqik8PFzjx4/Pt53z23+zZs00YcIE1alTR+XLl5fFYsn2VGCLxaK+ffsal1rv1lth\ns7LYzHqMGwCUkEOHDqlMmTJ68MEHJd0czDdt2jTFxMTIYrHoxRdfNB5ABtwt3nnnHX3wwQdydnbW\n5s2b7cbY3I24pALgrvfDDz/ozTfflKurq9zd3ZWcnKzr16/LYrGofv36GjBggKOrCBRYRESE4uPj\ndfLkSVksFj377LN3fdiQCBwA7gFNmjRR27ZttX//fiUnJ6t8+fJq0KCBOnbsqH79+tndGgrc6U6c\nOKFTp06pWrVq6tq1q1555RVHV6lYcEkFAACYjkGjAADAdAQOAABgOgIHAAAwHYEDAACYjsABAABM\nR+AAAACm+/+TZNRuwmiT7gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6b89e51f50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "paper.hyper_figure_label_printer('hr_missense_by_livermet_pfs')\n",
    "sb.boxplot(data=livermet['grp_coefs'], x='exp(beta)', y='group', fliersize=0, whis=[2.5, 97.5])\n",
    "_ = plt.xlim([0, 2])\n",
    "_ = plt.vlines(1, -10, 10, linestyles='--')\n",
    "#_ = plt.title('Hazard associated with log(missense SNV count / MB) \\n by liver-met status')\n",
    "_ = plt.ylabel('')\n",
    "_ = plt.xlabel('HR for {}'.format(cohort.hazard_plot_name))\n",
    "_ = plt.yticks([0, 1], ['No Liver Metastasis','With Liver Metastasis'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{{{hr_missense_pfs_livermet_0:n=16, HR=0.73, 95% CI (0.50, 1.02)}}}\n",
      "{{{hr_missense_pfs_livermet_1:n=9, HR=0.96, 95% CI (0.66, 1.37)}}}\n"
     ]
    }
   ],
   "source": [
    "livermet['grp_coefs']['group_label'] = livermet['grp_coefs'].group.apply(\n",
    "    lambda x: 'hr_missense_pfs_livermet_{}'.format(x)\n",
    ")\n",
    "for (name, livermet_class), group in livermet['grp_coefs'].groupby(['group_label','group']):\n",
    "    group_n = len(df.loc[df.liver_mets==livermet_class, 'liver_mets'])\n",
    "    paper.hyper_label_printer(formatter=paper.hr_posterior_formatter,\n",
    "                              label=name,\n",
    "                              n=group_n,\n",
    "                              series=group['exp(beta)'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{{{pval_missense_pfs_interaction_livermet:0.14}}}\n"
     ]
    }
   ],
   "source": [
    "## bayesian p-value for interaction\n",
    "comparison2 = pd.pivot_table(livermet['grp_coefs'],\n",
    "                  index = ['iter', 'model_cohort', 'variable'],\n",
    "                  values = 'value', columns = 'group_label')\n",
    "pval=1-(comparison2.eval('hr_missense_pfs_livermet_0 < hr_missense_pfs_livermet_1').mean())\n",
    "print('{{{%s:%s}}}' % ('pval_missense_pfs_interaction_livermet',\n",
    "                       paper.float_str(pval)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  },
  "name": "TCR Clonality.ipynb",
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true
  },
  "toc_position": {
   "height": "877px",
   "left": "0px",
   "right": "2073px",
   "top": "106px",
   "width": "538px"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
